Nonequilibrium steady states in sheared binary fluids 
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We simulate by lattice Boltzmann the steady shearing of a binary fluid mixture undergoing phase 
separation with full hydrodynamics in two dimensions. Contrary to some theoretical scenarios, a 
dynamical steady state is attained with finite domain lengths La:,y in the directions {x,y) of velocity 



and velocity gradient. Apparent scaling exponents are estimated as 
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We discuss the relative roles of diffusivity and hydrodynamics in attaining steady state. 
PACS numbers: 47.11.+j 
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Systems that are not in thermal equilibrium play a 
central role in modern statistical physics, and arise in ar- 
eas ranging from soap manufacture to subcellular biology 
0. Such systems include two important classes: those 
that are evolving towards Boltzmann equilibrium (e.g., 
by phase separation following a temperature quench), 
and those that arc maintained in nonequilibrium by con- 
tinuous driving (such as a shear flow). Of fundamen- 
tal interest, and surprising physical subtlety, are systems 
combining both features — such as a binary fluid under- 
going phase separation in the presence of shear. Here it 
is not known ^ whether coarsening continues indefi- 
nitely, as it does without shear, or whether a steady state 
is reached, in which the characteristic length scales Lx,y,z 
of the fluid domain structure attain finite 7-dcpcndcnt 
values at late times. (We define the mean velocity as 
Ux = iy so that X, ?/, z are velocity, velocity gradient and 
vorticity directions respectively; 7 is the shear rate.) 

Experimentally, saturating length scales are reportedly 
reached after a period of anisotropic domain growth j|, ||. 
However, the extreme elongation of domains along the 
flow direction means that, even in experiments, finite 
size effects could play an essential role in such saturation 
Theories in which the velocity does not fluctuate, 
but does advect the diffusive fluctuations of the concen- 
tration field, predict instead indefinite coarsening, with 
length scales Ly z scaling as 7-independent powers of the 



time t since quench, and (typically) Lx ^ jtLy 



real fluids, however, the velocity fluctuates strongly in 
nonlinear response to the advected concentration fleld, 
and hydrodynamic scaling arguments, balancing either 
interfacial and viscous or interfacial and inertial forces. 



predict saturation (e.g., L ~ 7 ^ or L 
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Given these experimental and theoretical differences of 
opinion, computer simulations of sheared binary fluids, 
with full hydrodynamics, are of major interest. 

The aforementioned scaling arguments cannot really 
distinguish one Cartesian direction from another, but 
even in theories that can do so, a two dimensional (2D) 
representation, suppressing z, is expected to capture the 
main physics Q. (Without shear, subtle non-scaling 
effects arise in 2D from the formation of disconnected 



droplets but shear seems to suppress these |^.) Per- 
forming simulations in 2D is therefore a fair compromise, 
especially given the extreme computational demands of 
the full 3D problem |, But, apart from ||, most 
numerical studies of binary fluids under ste ady shear, 
even in 2D, neglect hydrodynamics altogether III4 [l3| . 
Among fully hydrodynamic simulations (e.g., ]9|, |ll[), 
only Wagner and Ycomans |^ make a strong case for dy- 
namical steady states. In some cases these authors found 
complete remixing of the fluids {L 0); in the remain- 
der, finite size effects could not be excluded. (To do so 
requires L^^y ^x,y for a A:^ x Ay simulation box.) The 
existence of nonequilibrium steady states, with finite ix,y 
in an infinite system, therefore remains an open question. 

In this letter wc extend the hydrodynamic lattice 
Boltzmann (LB) studies of Refs.|^, ^ to much larger 
systems, which we then study over several decades of non- 
dimensionalized shear rate. Unlike previous authors, we 
are able to give clear evidence of true dynamical steady 
states, uncontaminated by finite size effects or other arti- 
facts. (Finite size effects typically result in quasi-laminar 
stripe domains which connect with themselves after one 
or more circuits of the periodic boundary conditions 
1^, |l^.) We then combine datasets using a quantitative 
scaling methodology developed for the unsheared prob- 
lem in |l5| ; this allows scaling exponents to be estimated 
using combined multi-decade fits. By this method we 



In find apparent scaling exponents L.^. 
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sustained over six decades of shear rate. 



Our basic LB algorithm for binary fluids is essentially 
as reported in (see also jl6j) on a D2Q9 lattice. Ad- 
ditionally we exploit recent algorithmic advances 
that overcome the intrinsic fluid velocity limit of LB by 
using blockwise translating lattice slabs connected by 
multiple Lees- Edwards boundary conditions (De- 
tails of our boundary conditions, with validation data, 
appear in p^.) One technical problem that remains 
within our LB scheme concerns the role of order param- 
eter diffusivity. In the hydrodynamic coarsening regimes 
of main interest (late times, modest shear rates) this dif- 
fusivity should always maintain local equilibrium across 
fluid interfaces, but never transport significant material 
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Figure 1: Snapshots of the steady-state order-parameter for 
R020 with 7 — 0.0005 at dimensionless times t'y — 60 and 255 
(upper and lower). 



across the interior of domains ||15| . Under shear when 
domains are extremely anisotropic, compromise becomes 
inevitable. We discuss below the implications of this for 
the interpretation of our apparent scaling exponents. 

The governing equations that our LB scheme ap- 
proximates are the Cahn-Hilliard equation for the com- 
positional order parameter ip, and the incompressible 
(VqUq = 0) Navier-Stokes equation for the velocity Ua 
in an isothermal fluid of unit mass density: 

{dtUa + UpV pUa) + \7aP — l^^^Ua - (f>Vafi = 0(1) 
dt(p + Va {(fUc - MVafi) = (2) 

Here p is pressure (related in LB to density fluctua- 
tions, which remain small, by an ideal gas equation of 
15 ); I' is the kinematic viscosity; M is the order- 



state 

parameter mobility, and ^ = Bip [ip^ — l) — kV"^^, the 
chemical potential. B and k are positive constants; the 
interfacial tension is a = (SkB/O)^/^ and the interfacial 
width is ^0 = (2k/_B)^/^ LB control parameters are 
B, K, M and v alongside the steady shear rate 7. 

Blockwise sheared boundary conditions are imposed 
[0 such that " Vua; dy = Aj,7. A fluctuating local 
velocity field can then arise by nonlincarity, as in experi- 
ments ll^. (We neglect thermal fluctuations in our fluid, 
as appropriate for dynamics near a zero-temperature 
fixed point [^.) Under shear we define length scales 
L,j. y using a gradient statistic for ip that measures the 
mean distance between interfaces lying across to the cho- 
sen direction B. We also define L\\ ^. with Lii > L±^ by 



reference to appropriate principal axes. (Ref. finds 
universality, in a related but nonhydrodynamic system, 
when lengths are scaled with but not with Lx,y.) 

The physics of binary fluid demixing, with no shear 
and low enough diffusivity A/, can be nondimcnsional- 
ized via a single length scale Lq = v j a and time scale 
To = ja^. These are, up to dimensionless prcfactors 
|]l5| , the domain size and time after quench at which 
an interfacial/ viscous balance in the coarsening dynam- 
ics (viscous hydrodynamic regime, VH) crosses over to an 
intcrfacial/inertial balance (inertia! hydrodynamic, IH). 
By scaling the domain length L and times t by Lq and Tqi 
multiple datasets were shown to merge, giving a universal 
crossover between VH and IH regimes Q . 

Accordingly, in our search for nonequilibrium steady 
states, we nondimensionalize the shear rate as ^Tq and 
domain sizes as y || j^/Lq. One can then expect plots 
of length against strain rate (or its inverse, as used below) 
to show data collapse whenever diffusivity is small. One 
possibility ||^ is that all such plots might exhibit the same 
crossover from VH to IH behaviour as would be found by 
substituting t = 7"^ in the universal scaling plot of L/ Lq 
against t/T^. This would give Lx^y^\\^±_ ^ L with L/ ^ 
(7To)~^ at large shear rates and L/Lq ~ (t'Tq)"^/'^ at 
small ones. Alternatively, some single power law could 
persist at all ^Tq 0|; and/or there could be different 
exponents in the different directions; or there might be 
no steady state at all ||] . 

All simulations reported here were done for fully sym- 
metric quenches on a A^, x Ay = 1024 x 512 lattice, with 
up to i = 6 X 10^ updates. Parameters ^Ojf, M, i/ (Ta- 
ble |) and 7 were chosen, following ||T^, so that: inter- 
faces are wide enough to be resolved (restrictions on ^0); 
fluid flow is slow enough for advected interfaces to be in 
local equilibrium (restrictions on <j and 7); the diffusiv- 
ity is low enough not to contaminate steady-state length 
scales, e.g., detectable as a strong residual M dependence 
(restriction on M). Also, these steady-state lengths must 
be sufficiently small to avoid finite size effects (quantified 
below) and the code must run stably for long enough to 
achieve steady state, typically t > 10^ updates. 

Acceptable shear rates were found to be 1.25 x 10~^ < 
7 < 2 X 10~^ (in lattice units). Higher values gave inac- 
curacies as listed above; lower values gave unacceptably 
long run times. As in the unsheared case |l^ judicious 
combinations of f , M and v allow systems spanning 
several decades in L/ and 7T0 to be accurately studied, 
by exploiting LB's ability to vary Lq and Tq alongside 7. 

Fig. 1^ shows two snapshots of the order-parameter field 
for R020 with 7 = 5 x lO"'* after a steady state had 
been reached. The snapshots are at dimensionless times 
t7 = 60 and t7 = 255, for the upper and lower plots 
respectively. Fig. ^ shows unsealed time-series for L^ 
and Ly from a representative set of simulations with 7 = 
5 X 10~^. Both figures show decisive evidence of length- 
scale saturation in a regime that seems safely clear of 
any finite size effects. A number of tests were performed 
in which all run parameters were held constant except 
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Name 


V 


M 


'-^theory 




u 


To 


R028 


1.41 


0.05 


0.063 


0.055 


36.1 


927 


R022 


0.5 


0.25 


0.047 


0.042 


5.95 


70.9 


R029 


0.2 


0.15 


0.047 


0.042 


0.952 


4.54 


R020 


0.025 


2 


0.0047 


0.0042 


0.149 


0.886 


R030 


0.00625 


1.25 


0.0047 


0.0042 


0.00930 


0.0138 


R019 


0.0014 


4 


0.0024 


0.0021 


0.000933 


0.000622 


R032 


0.0005 


5 


0.00094 


0.00083 


0.000301 


0.000181 



Table I: Parameter sets used in simulations, along with Lo 
and Tq. In all cases = 1-13. (See |l5j for discussion on the 
difference between the theoretical and measured values of cr; 
the measured ones are used to determine Lo,To). 



o 

+ 



CO 

o 



o 

+ 







mm*/ 





"1 I I r 

1e+01 1e+03 1e+05 1e+07 




-A 




„-,..',"'"-,..s"..'.'rV'''V,-"'"'"-'-,.''.,»'R028 - 

R02G - 

R030 -■ 

R019 - 



t (X 10^) 



70 
60 
50 
40 
30 
20 
10 





R028 
R020 
R030 
R019 



/ (X 10^) 



Figure 2: Plots of Lx{t) and Ly{t) in lattice units, for various 
runs with 7 = 5 X lO"**. 



the lattice dimensions which were changed in the ranges 



A, 



512 to 2048 and A.„ = 256 to 1024. From the results 



of these tests we conclude that finite-size effects are fully 
under control when L^^y < A^^y/i, a criterion extensively 
benchmarked in unsheared systems ||l^. At the same 
time, ^ 30 in lattice units, well clear of discretisation 
artifacts. However, the thin fluid threads visible in Fig. |l| 
mean that residual diffusion cannot entirely be ruled out; 
we return to this point below. 
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Figure 3: Dimensionless scaling plots of lengths vs shear rate. 
Error bars as shown give the 95% confidence limits (dashed 
lines) on the fitted linear regression (solid line). Solid symbols: 
7 = 5 X 10~*. Open symbols have 7 one of, from left to right, 
2 X 10-^ 10-^ 2.5 X 10"* or 1.25 x 10"*. The fit used only 
the data for 7 = 5 x 10"*. Points with two symbols at the 
same l/7ro denote variations in the diffusivity M. 



Fig. H shows dimensionless scaling plots of steady-state 
length scales L^.y against shear rate, for a series of runs 
in which 7 = 5 x 10~* with variable Lq,Tq (solid sym- 
bols). Lx,y were obtained as the temporal means of the 
time-series data Lx{t) and Ly{t), after discarding data 
for which t < 10^. Uncertainties in L^^y (also i||,±) were 
found using 999 bootstrap replicas of the time-series data. 
These uncertainties, which provide the relative error bars 
indicated in the scaling plots, were then used to weight 
the individual data-points in a linear least-squares regres- 
sion from which estimates were obtained of the intercept 
and gradient of the straight lines shown on the log-log 
plots. Fitted exponents for and are —0.678 ±0.042 
and -0.759 ± 0.029; and for L|| and Lj_, -0.678 ± 0.039 
and -0.756 ±0.030 (data not shown) ||l|. The 95% con- 
fidence level admits the appealing ansatz of fractional 



power laws 
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However, eaution is warranted in presenting these ap- 
parently clean scaling laws. Firstly, were the dynamic 
length scales to be found by substituting t ~ 'Y~-^ in the 
shear- free coarsening plot as suggested above ||^ , then a 
slow crossover from VH (L ^ 7^^) to IH {L ~ 7^^^'^) 
would affect this entire range of 'jTo [|l5|. Fitting L^^y 
data to simple power laws might therefore be misleading. 
Secondly, we also show in Fig. js] two datasets found by 
varying 7 with other parameters fixed. For both and 
Ly, if taken in isolation these sets would suggest smaller 
slopes than seen for the main plot. Such deviations were 
found previously in the unsheared case jist , and argued 
to be a signature of residual diffusion, with each dataset 
asymptoting onto the global trend line from above left. 

The upward curvature of these two datasets, and the 
fact that they appear to asymptote onto (rather than 
cross) the best fit line found for 7 = 5 x 10^"', offers some 
reassurance, but no guarantee, that the latter is uncon- 
taminated by residual diffusion at the domain scale. If so, 
our stationary states stem not from diffusion (although 
that would itself be interesting) but from the hydrody- 
namic balance between the stretching, breaking and coa- 
lescence of domains. Further evidence for this comes from 
study of the evolving ip field, in which all of these effects 
are visible but (with our parameters) large-scale diffusion 
is not. Also, all the simulations reported here quench at 
t = from a noisy but uniform state. The system then 
passes through, and apparently leaves, a diffusive regime 
prior to the hydrodynamic one that gets cut off by the 
presence of shear. As a final check, we show directly the 
effect of varying M in two runs shown in Fig. (where 



more than one symbol occurs at the same 7T0). The 
lower data-points were found by reducing M by a fac- 
tor 2 or more from nearby runs. The resulting shifts are 
modest, although not entirely negligible — particularly 
towards the bottom left of the plots. Thus it remains pos- 
sible that further reduction in M (not practical numeri- 
cally at present) could reveal a significant kink on these 
plots, as might be expected near a VH to IH crossover. 

In conclusion, although the apparent scaling exponents 
reported above are interesting and merit both theoreti- 
cal investigation and experimental tests, the main signif- 
icance of our work is in the unambiguous demonstration 
of nonequilibrium steady states in sheared binary fluids. 
Since theories that neglect velocity fluctuations do not 
predict such states hydrodynamics appears to play 
an essential role. (This remains true even if diffusion is 
not negligible as considered above.) A key question is 
whether such steady states persist in three dimensions. 
Although the physics of stretching, breakup and coales- 
cence is captured in 2D, in 3D there can remain tubular 
connections between domains in the vorticity direction; 
these could remain relatively unaffected by shear. This 
might leave open a route to continuous coarsening that 
is topologically absent in two dimensions. Wc hope to 
address the 3D case in future simulations. 
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